input data
source(file = '../inputProcessing/date_input.r')
rep <- 1e4
rep_sim <- 1e3
for (date_num in 1:length(dates_input)){
# date_num = 1
date_week_finishing <- as.Date(dates_input[date_num],format = '%Y-%m-%d')
output <- list()
output2 <- list()
for (Mdata in c('A','G')){
# Mdata = 'A'
inputs<- readRDS(file=paste0(
'../../Mobility.2.0.Save/Saved_file/inputProcessing/Rdata/',
date_week_finishing,'_inputs_from_',Mdata,'.rds'))
for (si in 1:2){
# si=1
########################################
####inputs
D <- inputs$D
M <- inputs$M
Ot <- inputs$Ot[[si]]$Ot
W <- inputs$W[[si]]$W
H <- inputs$H
SI <- inputs$SI[[si]]$SI
delta_id <- inputs$delta_id
#
country <- names(D)[-1]
mD <- as.matrix(D[,-1])
N_geo <- ncol(D)-1
N_days <- nrow(D)
N_week <- nrow(D)/7
# prerp MCMC
# epi
R0 <- 5
# for risk
b <- 0
o <- 1
theta0 <- c(rep(R0,N_geo),
rep(b,N_geo),
rep(o,N_geo))
prior_theta <- matrix(c(rep(c(0,5),N_geo),
rep(c(-1e2,1e2),N_geo),
rep(c(0,1e5),N_geo)),
length(theta0),2, byrow=TRUE)
# parameter names
f0 <- function(x) paste0('R0_',x)
f1 <- function(x) paste0('beta_',x)
f2 <- function(x) paste0('Over_',x)
n_t<- c(sapply(country,f0), sapply(country,f1), sapply(country,f2))
# sd dev for proposal
sigma <- rep(1e-1,length(theta0))
########################################
# useful functions
# useful functions
sapply(paste0('Rscript/std/',(list.files('Rscript/std/'))),FUN = source)
########################################
#run MCMC
#check
# res <- MCMC_iter(iter = rep, theta0 = theta0, s = sigma)
res <- MCMC_full(iter = rep,
theta0 = theta0,
s = sigma,
repli_adapt = 10,
within_iter = rep/10)
########################################
## check convergence
print('########################################')
print(paste0('check convergence',date_week_finishing,' - ',Mdata,' - ',si))
Acc <- colSums(diff(res$theta)!=0)/(rep-1)
if(rep >1e2) f <- round(seq(1,rep,length.out = 1e2))
plot(res$logL[,1],main=paste0('DIC ',res$DIC[1],', P ',res$DIC[2]))
layout(matrix(1:4,2,2))
for (i in 1:length(theta0)){
plot(res$theta[,i],
main = paste0(n_t[i],' - ',round(Acc[i]*100)))#, ylim = prior_theta[i,])
}
# apply(res$theta,2,quantile,c(.5,.025,.975))
ncountry <- country
ncountry[which(ncountry %in% 'United_Kingdom')] <- 'UK'
ncountry[which(ncountry %in% 'United_States_of_America')] <- 'USA'
R0s <- apply(res$theta[,1:N_geo],2,quantile,c(.5,.025,.975))
betas <- apply(res$theta[,N_geo+(1:N_geo)],2,quantile,c(.5,.025,.975))
over <- apply(res$theta[,2*N_geo+(1:N_geo)],2,quantile,c(.5,.025,.975))
## margin for side 2 is 7 lines in siz
layout(1)
par(mar = c(7,4,4,2) + 0.1) ## default is c(5,4,4,2) + 0.1
errbar(1:N_geo,R0s[1,],R0s[2,],R0s[3,],
xlab = '', ylab = TeX('R_0'), bty = 'n',xaxt = "n",ylim=c(0,10))
# xlab = '', ylab = 'R',ylim = c(0,3), bty = 'n',xaxt = "n")
lines(c(1,N_geo),rep(1,2), col = 'red')
axis(1, at=1:N_geo, labels=ncountry,las=2)
####
errbar(1:N_geo,betas[1,],betas[2,],betas[3,],
xlab = '', ylab = TeX('$\\beta$'), bty = 'n',xaxt = "n")
# xlab = '', ylab = 'R',ylim = c(0,3), bty = 'n',xaxt = "n")
# lines(c(1,N_geo),rep(1,2), col = 'red')
axis(1, at=1:N_geo, labels=ncountry,las=2)
output[[si]] <- list(resMCMC = res,
R0s = R0s,
betas = betas,
over = over
)
####################################
## Rt daily and Rt_D2
### observed mobility
res_m <- output[[si]]
rep <- nrow(res_m$resMCMC$theta)
Rt_daily <- array(NA,dim = c(N_days,N_geo,rep))
Rt_D2 <- Rt_daily
m_eff <- Rt_daily
f_m_eff <- function(){
1 + 1/B * log( inputs$H %*% exp( - B * (1-inputs$M) ) )
}
for (j in 1:rep){
R_daily <- Rt_fun(theta = res_m$resMCMC$theta[j,], x = inputs$M )
Rt_daily[,,j] <- R_daily
Rt_D2[,,j] <- inputs$H %*% R_daily
B <- rep(1,nrow(inputs$M)) %*% t(res_m$resMCMC$theta[j,(2-1)*N_geo+ (1:N_geo)])
m_eff[,,j] <- f_m_eff() # 1 + 1/B * log( inputs$H %*% exp( - B * (1-inputs$M) ) )
}
###################################
### assumed mobility
n_d <- 1e3
Rt_assumed_mob <- array(NA,dim = c(n_d,N_geo,rep))
x <- matrix(seq(0,1,length.out = n_d),n_d,N_geo,byrow = FALSE)
for (j in 1:rep){
R_daily <- Rt_fun(theta = res_m$resMCMC$theta[j,], x = 1-x )
Rt_assumed_mob[,,j] <- R_daily
}
#####################################
## summaries
### Summary for Rt daily and Rt_D2
# format outputs
temp <- D
temp[,-1] <- NA
results_full_Rt_daily <- list(median_R = temp,
low_R = temp,
up_R = temp,
M = temp)
results_full_Rt_D2 <- list(median_R = temp,
low_R = temp,
up_R = temp,
M_D = temp)
results_meff <- list(median_meff = temp,
low_meff = temp,
up_meff = temp)
for (i in 1:N_geo){
temp <- apply(Rt_daily[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
results_full_Rt_daily$median_R[,i+1] <- temp[1,]
results_full_Rt_daily$low_R[,i+1] <- temp[2,]
results_full_Rt_daily$up_R[,i+1] <- temp[3,]
results_full_Rt_daily$M[,i+1] <- M[,i]
temp <- apply(Rt_D2[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
results_full_Rt_D2$median_R[,i+1] <- temp[1,]
results_full_Rt_D2$low_R[,i+1] <- temp[2,]
results_full_Rt_D2$up_R[,i+1] <- temp[3,]
# results_full_Rt_D2$M_D[,i+1] <- M_D[,i]
temp <- apply(m_eff[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
results_meff$median_meff[,i+1] <- temp[1,]
results_meff$low_meff[,i+1] <- temp[2,]
results_meff$up_meff[,i+1] <- temp[3,]
}
#############################################
### Summary for Rt assumed mobility
# format outputs
temp <- data.frame(matrix(NA,nrow = n_d, ncol = 1+N_geo))
temp[,1] <- x[,1]
names(temp) <- c('mobility', country)
results_full_Rt_assumed_mob <- list(median_R = temp,
low_R = temp,
up_R = temp)
Thresholds <- output[[1]]$R0s
f <- c()
for (i in 1:N_geo){
temp <- apply(Rt_assumed_mob[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
results_full_Rt_assumed_mob$median_R[,i+1] <- temp[1,]
results_full_Rt_assumed_mob$low_R[,i+1] <- temp[2,]
results_full_Rt_assumed_mob$up_R[,i+1] <- temp[3,]
f[1] <- which(results_full_Rt_assumed_mob$median_R[,i+1] < 1)[1]
f[2] <- which(results_full_Rt_assumed_mob$low_R[,i+1] < 1)[1]
f[3] <- which(results_full_Rt_assumed_mob$up_R[,i+1] < 1)[1]
for(j in 1:3){
if (is.na(f[j])){
Thresholds[j,i] <- NA
}else{
Thresholds[j,i] <-mean(results_full_Rt_assumed_mob$median_R$mobility[f[j]:(f[j]+1)])
}
}
}
output[[si]]$Thresholds <- Thresholds
#############################################
# forecast and predictions
## short-term
n_pred <- 7
# choose future mobility
M_pred <- rbind(M,matrix(M[nrow(M),],n_pred,N_geo,byrow = TRUE))
# useful functions
sapply(paste0('Rscript/proj_f/',(list.files('Rscript/proj_f'))),FUN = source)
res_projection_ST <- mobility_prediction_over2(n_pred = n_pred, rep_sim = rep_sim)
## summary
temp <- data.frame(dates = seq(1,n_pred) + D$dates[nrow(D)],
matrix(NA,n_pred,N_geo))
names(temp) <- c('dates',country)
summary_proj <- list(median = temp,
low = temp,
up = temp)
for (i in 1:N_geo){
temp <- apply(res_projection_ST$D_pred[,i,],1,quantile,c(.5,.025,.975))
summary_proj$median[,i+1] <- temp[1,]
summary_proj$low[,i+1] <- temp[2,]
summary_proj$up[,i+1] <- temp[3,]
}
output2[[si]] <- list(results_full_Rt_daily = results_full_Rt_daily,
results_full_Rt_D2 = results_full_Rt_D2,
results_meff = results_meff,
results_full_Rt_assumed_mob = results_full_Rt_assumed_mob,
summary_proj_ST = summary_proj)
}
saveRDS(object = output,
file = paste0('../../Mobility.2.0.Save/Saved_file/FitModels/R0_beta_Over/Rdata/',
date_week_finishing,'_output_over_form_',Mdata,'.rds' ))
saveRDS(object = output2,
file = paste0('../../Mobility.2.0.Save/Saved_file/FitModels/R0_beta_Over/Rdata/',
date_week_finishing,'_output2_over_form_',Mdata,'.rds' ))
}
}
## [1] "########################################"
## [1] "check convergence2020-05-03 - A - 1"
## [1] "########################################"
## [1] "check convergence2020-05-03 - A - 2"
## [1] "########################################"
## [1] "check convergence2020-05-03 - G - 1"
## [1] "########################################"
## [1] "check convergence2020-05-03 - G - 2"
## [1] "########################################"
## [1] "check convergence2020-05-10 - A - 1"
## [1] "########################################"
## [1] "check convergence2020-05-10 - A - 2"
## [1] "########################################"
## [1] "check convergence2020-05-10 - G - 1"
## [1] "########################################"
## [1] "check convergence2020-05-10 - G - 2"
